home *** CD-ROM | disk | FTP | other *** search
/ Aminet 1 (Walnut Creek) / Aminet - June 1993 [Walnut Creek].iso / usenet / sources / volume2 / aplictns / matlab / src.5 < prev    next >
Internet Message Format  |  1988-11-02  |  56KB
  1. Path: xanth!nic.MR.NET!hal!cwjcc!mailrus!ulowell!page
  2. From: page@swan.ulowell.edu (Bob Page)
  3. Newsgroups: comp.sources.amiga
  4. Subject: v02i045:  matlab - matrix laboratory, Part05/11
  5. Message-ID: <10020@swan.ulowell.edu>
  6. Date: 2 Nov 88 21:44:05 GMT
  7. Organization: University of Lowell, Computer Science Dept.
  8. Lines: 1220
  9. Approved: page@swan.ulowell.edu
  10.  
  11. Submitted-by: strovink%galaxy-43@afit-ab.arpa (Mark A. Strovink)
  12. Posting-number: Volume 2, Issue 45
  13. Archive-name: applications/matlab/src.5
  14.  
  15. #    This is a shell archive.
  16. #    Remove everything above and including the cut line.
  17. #    Then run the rest of the file through sh.
  18. #----cut here-----cut here-----cut here-----cut here----#
  19. #!/bin/sh
  20. # shar:    Shell Archiver
  21. #    Run the following text with /bin/sh to create:
  22. #    src-5
  23. # This archive created: Wed Nov  2 16:20:35 1988
  24. cat << \SHAR_EOF > src-5
  25.             IF (N .LT. KP1) GO TO 100     
  26.             DO 90 J = KP1, N              
  27.                TR = AR(K,J)               
  28.                TI = AI(K,J)               
  29.                AR(K,J) = 0.0D0            
  30.                AI(K,J) = 0.0D0            
  31.                CALL WAXPY(K,TR,TI,AR(1,K),AI(1,K),1,AR(1,J),AI(1,J),1)          
  32.    90       CONTINUE   
  33.   100       CONTINUE   
  34.   110    CONTINUE      
  35. C                      
  36. C        FORM INVERSE(U)*INVERSE(L)       
  37. C                      
  38.          NM1 = N - 1   
  39.          IF (NM1 .LT. 1) GO TO 150        
  40.          DO 140 KB = 1, NM1               
  41.             K = N - KB                    
  42.             KP1 = K + 1                   
  43.             DO 120 I = KP1, N             
  44.                WORKR(I) = AR(I,K)         
  45.                WORKI(I) = AI(I,K)         
  46.                AR(I,K) = 0.0D0            
  47.                AI(I,K) = 0.0D0            
  48.   120       CONTINUE   
  49.             DO 130 J = KP1, N             
  50.                TR = WORKR(J)              
  51.                TI = WORKI(J)              
  52.                CALL WAXPY(N,TR,TI,AR(1,J),AI(1,J),1,AR(1,K),AI(1,K),1)          
  53.   130       CONTINUE   
  54.             L = IPVT(K)                   
  55.             IF (L .NE. K)                 
  56.      *         CALL WSWAP(N,AR(1,K),AI(1,K),1,AR(1,L),AI(1,L),1)                
  57.   140    CONTINUE      
  58.   150    CONTINUE      
  59.   160 CONTINUE         
  60.       RETURN           
  61.       END              
  62.       SUBROUTINE WPOFA(AR,AI,LDA,N,INFO)  
  63.       DOUBLE PRECISION AR(LDA,1),AI(LDA,1)                   
  64.       DOUBLE PRECISION S,TR,TI,WDOTCR,WDOTCI                 
  65.       DO 30 J = 1, N   
  66.          INFO = J      
  67.          S = 0.0D0     
  68.          JM1 = J-1     
  69.          IF (JM1 .LT. 1) GO TO 20         
  70.          DO 10 K = 1, JM1                 
  71.             TR = AR(K,J)-WDOTCR(K-1,AR(1,K),AI(1,K),1,AR(1,J),AI(1,J),1)        
  72.             TI = AI(K,J)-WDOTCI(K-1,AR(1,K),AI(1,K),1,AR(1,J),AI(1,J),1)        
  73.             CALL WDIV(TR,TI,AR(K,K),AI(K,K),TR,TI)           
  74.             AR(K,J) = TR                  
  75.             AI(K,J) = TI                  
  76.             S = S + TR*TR + TI*TI         
  77.    10    CONTINUE      
  78.    20    CONTINUE      
  79.          S = AR(J,J) - S                  
  80.          IF (S.LE.0.0D0 .OR. AI(J,J).NE.0.0D0) GO TO 40      
  81.          AR(J,J) = DSQRT(S)               
  82.    30 CONTINUE         
  83.       INFO = 0         
  84.    40 RETURN           
  85.       END              
  86.       SUBROUTINE RREF(AR,AI,LDA,M,N,EPS)  
  87.       DOUBLE PRECISION AR(LDA,1),AI(LDA,1),EPS,TOL,TR,TI,WASUM                  
  88.       TOL = 0.0D0      
  89.       DO 10 J = 1, N   
  90.          TOL = DMAX1(TOL,WASUM(M,AR(1,J),AI(1,J),1))         
  91.    10 CONTINUE         
  92.       TOL = EPS*DFLOAT(2*MAX0(M,N))*TOL   
  93.       K = 1            
  94.       L = 1            
  95.    20 IF (K.GT.M .OR. L.GT.N) RETURN      
  96.       I = IWAMAX(M-K+1,AR(K,L),AI(K,L),1) + K-1              
  97.       IF (DABS(AR(I,L))+DABS(AI(I,L)) .GT. TOL) GO TO 30     
  98.          CALL WSET(M-K+1,0.0D0,0.0D0,AR(K,L),AI(K,L),1)      
  99.          L = L+1       
  100.          GO TO 20      
  101.    30 CALL WSWAP(N-L+1,AR(I,L),AI(I,L),LDA,AR(K,L),AI(K,L),LDA)                 
  102.       CALL WDIV(1.0D0,0.0D0,AR(K,L),AI(K,L),TR,TI)           
  103.       CALL WSCAL(N-L+1,TR,TI,AR(K,L),AI(K,L),LDA)            
  104.       AR(K,L) = 1.0D0                     
  105.       AI(K,L) = 0.0D0                     
  106.       DO 40 I = 1, M   
  107.          TR = -AR(I,L)                    
  108.          TI = -AI(I,L)                    
  109.          IF (I .NE. K) CALL WAXPY(N-L+1,TR,TI,               
  110.      $                 AR(K,L),AI(K,L),LDA,AR(I,L),AI(I,L),LDA)                 
  111.    40 CONTINUE         
  112.       K = K+1          
  113.       L = L+1          
  114.       GO TO 20         
  115.       END              
  116.       SUBROUTINE HILBER(A,LDA,N)          
  117.       DOUBLE PRECISION A(LDA,N)           
  118. C     GENERATE INVERSE HILBERT MATRIX     
  119.       DOUBLE PRECISION P,R                
  120.       P = DFLOAT(N)    
  121.       DO 20 I = 1, N   
  122.         IF (I.NE.1) P = (DFLOAT(N-I+1)*P*DFLOAT(N+I-1))/DFLOAT(I-1)**2          
  123.         R = P*P        
  124.         A(I,I) = R/DFLOAT(2*I-1)          
  125.         IF (I.EQ.N) GO TO 20              
  126.         IP1 = I+1      
  127.         DO 10 J = IP1, N                  
  128.           R = -(DFLOAT(N-J+1)*R*(N+J-1))/DFLOAT(J-1)**2      
  129.           A(I,J) = R/DFLOAT(I+J-1)        
  130.           A(J,I) = A(I,J)                 
  131.    10   CONTINUE       
  132.    20 CONTINUE         
  133.       RETURN           
  134.       END              
  135.       SUBROUTINE HTRIDI(NM,N,AR,AI,D,E,E2,TAU)               
  136. C                      
  137.       INTEGER I,J,K,L,N,II,NM,JP1         
  138.       DOUBLE PRECISION AR(NM,N),AI(NM,N),D(N),E(N),E2(N),TAU(2,N)               
  139.       DOUBLE PRECISION F,G,H,FI,GI,HH,SI,SCALE               
  140.       DOUBLE PRECISION FLOP,PYTHAG        
  141. C                      
  142. C     THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF                 
  143. C     THE ALGOL PROCEDURE TRED1, NUM. MATH. 11, 181-195(1968)                   
  144. C     BY MARTIN, REINSCH, AND WILKINSON.  
  145. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).           
  146. C                      
  147. C     THIS SUBROUTINE REDUCES A COMPLEX HERMITIAN MATRIX     
  148. C     TO A REAL SYMMETRIC TRIDIAGONAL MATRIX USING           
  149. C     UNITARY SIMILARITY TRANSFORMATIONS.                    
  150. C                      
  151. C     ON INPUT.        
  152. C                      
  153. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL                 
  154. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM                  
  155. C          DIMENSION STATEMENT.           
  156. C                      
  157. C        N IS THE ORDER OF THE MATRIX.    
  158. C                      
  159. C        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,     
  160. C          RESPECTIVELY, OF THE COMPLEX HERMITIAN INPUT MATRIX.                 
  161. C          ONLY THE LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED.              
  162. C                      
  163. C     ON OUTPUT.       
  164. C                      
  165. C        AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY TRANS-                 
  166. C          FORMATIONS USED IN THE REDUCTION IN THEIR FULL LOWER                 
  167. C          TRIANGLES.  THEIR STRICT UPPER TRIANGLES AND THE  
  168. C          DIAGONAL OF AR ARE UNALTERED.  
  169. C                      
  170. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE THE TRIDIAGONAL MATRIX.        
  171. C                      
  172. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL                 
  173. C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO.              
  174. C                      
  175. C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.            
  176. C          E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.                
  177. C                      
  178. C        TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS.            
  179. C                      
  180. C     MODIFIED TO GET RID OF ALL COMPLEX ARITHMETIC, C. MOLER, 6/27/79.         
  181. C                      
  182. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,                
  183. C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY                 
  184. C                      
  185. C     ------------------------------------------------------------------        
  186. C                      
  187.       TAU(1,N) = 1.0D0                    
  188.       TAU(2,N) = 0.0D0                    
  189. C                      
  190.       DO 100 I = 1, N                     
  191.   100 D(I) = AR(I,I)   
  192. C     .......... FOR I=N STEP -1 UNTIL 1 DO -- ..........    
  193.       DO 300 II = 1, N                    
  194.          I = N + 1 - II                   
  195.          L = I - 1     
  196.          H = 0.0D0     
  197.          SCALE = 0.0D0                    
  198.          IF (L .LT. 1) GO TO 130          
  199. C     .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) ..........               
  200.          DO 120 K = 1, L                  
  201.   120    SCALE = FLOP(SCALE + DABS(AR(I,K)) + DABS(AI(I,K))) 
  202. C                      
  203.          IF (SCALE .NE. 0.0D0) GO TO 140  
  204.          TAU(1,L) = 1.0D0                 
  205.          TAU(2,L) = 0.0D0                 
  206.   130    E(I) = 0.0D0                     
  207.          E2(I) = 0.0D0                    
  208.          GO TO 290     
  209. C                      
  210.   140    DO 150 K = 1, L                  
  211.             AR(I,K) = FLOP(AR(I,K)/SCALE)                    
  212.             AI(I,K) = FLOP(AI(I,K)/SCALE)                    
  213.             H = FLOP(H + AR(I,K)*AR(I,K) + AI(I,K)*AI(I,K))  
  214.   150    CONTINUE      
  215. C                      
  216.          E2(I) = FLOP(SCALE*SCALE*H)      
  217.          G = FLOP(DSQRT(H))               
  218.          E(I) = FLOP(SCALE*G)             
  219.          F = PYTHAG(AR(I,L),AI(I,L))      
  220. C     .......... FORM NEXT DIAGONAL ELEMENT OF MATRIX T ..........              
  221.          IF (F .EQ. 0.0D0) GO TO 160      
  222.          TAU(1,L) = FLOP((AI(I,L)*TAU(2,I) - AR(I,L)*TAU(1,I))/F)               
  223.          SI = FLOP((AR(I,L)*TAU(2,I) + AI(I,L)*TAU(1,I))/F)  
  224.          H = FLOP(H + F*G)                
  225.          G = FLOP(1.0D0 + G/F)            
  226.          AR(I,L) = FLOP(G*AR(I,L))        
  227.          AI(I,L) = FLOP(G*AI(I,L))        
  228.          IF (L .EQ. 1) GO TO 270          
  229.          GO TO 170     
  230.   160    TAU(1,L) = -TAU(1,I)             
  231.          SI = TAU(2,I)                    
  232.          AR(I,L) = G   
  233.   170    F = 0.0D0     
  234. C                      
  235.          DO 240 J = 1, L                  
  236.             G = 0.0D0                     
  237.             GI = 0.0D0                    
  238. C     .......... FORM ELEMENT OF A*U ..........              
  239.             DO 180 K = 1, J               
  240.                G = FLOP(G + AR(J,K)*AR(I,K) + AI(J,K)*AI(I,K))                  
  241.                GI = FLOP(GI - AR(J,K)*AI(I,K) + AI(J,K)*AR(I,K))                
  242.   180       CONTINUE   
  243. C                      
  244.             JP1 = J + 1                   
  245.             IF (L .LT. JP1) GO TO 220     
  246. C                      
  247.             DO 200 K = JP1, L             
  248.                G = FLOP(G + AR(K,J)*AR(I,K) - AI(K,J)*AI(I,K))                  
  249.                GI = FLOP(GI - AR(K,J)*AI(I,K) - AI(K,J)*AR(I,K))                
  250.   200       CONTINUE   
  251. C     .......... FORM ELEMENT OF P ..........                
  252.   220       E(J) = FLOP(G/H)              
  253.             TAU(2,J) = FLOP(GI/H)         
  254.             F = FLOP(F + E(J)*AR(I,J) - TAU(2,J)*AI(I,J))    
  255.   240    CONTINUE      
  256. C                      
  257.          HH = FLOP(F/(H + H))             
  258. C     .......... FORM REDUCED A ..........                   
  259.          DO 260 J = 1, L                  
  260.             F = AR(I,J)                   
  261.             G = FLOP(E(J) - HH*F)         
  262.             E(J) = G   
  263.             FI = -AI(I,J)                 
  264.             GI = FLOP(TAU(2,J) - HH*FI)   
  265.             TAU(2,J) = -GI                
  266. C                      
  267.             DO 260 K = 1, J               
  268.                AR(J,K) = FLOP(AR(J,K) - F*E(K) - G*AR(I,K)   
  269.      X        + FI*TAU(2,K) + GI*AI(I,K)) 
  270.                AI(J,K) = FLOP(AI(J,K) - F*TAU(2,K) - G*AI(I,K)                  
  271.      X        - FI*E(K) - GI*AR(I,K))     
  272.   260    CONTINUE      
  273. C                      
  274.   270    DO 280 K = 1, L                  
  275.             AR(I,K) = FLOP(SCALE*AR(I,K))                    
  276.             AI(I,K) = FLOP(SCALE*AI(I,K))                    
  277.   280    CONTINUE      
  278. C                      
  279.          TAU(2,L) = -SI                   
  280.   290    HH = D(I)     
  281.          D(I) = AR(I,I)                   
  282.          AR(I,I) = HH                     
  283.          AI(I,I) = FLOP(SCALE*DSQRT(H))   
  284.   300 CONTINUE         
  285. C                      
  286.       RETURN           
  287.       END              
  288.       SUBROUTINE HTRIBK(NM,N,AR,AI,TAU,M,ZR,ZI)              
  289. C                      
  290.       INTEGER I,J,K,L,M,N,NM              
  291.       DOUBLE PRECISION AR(NM,N),AI(NM,N),TAU(2,N),ZR(NM,M),ZI(NM,M)             
  292.       DOUBLE PRECISION H,S,SI,FLOP        
  293. C                      
  294. C     THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF                 
  295. C     THE ALGOL PROCEDURE TRBAK1, NUM. MATH. 11, 181-195(1968)                  
  296. C     BY MARTIN, REINSCH, AND WILKINSON.  
  297. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).           
  298. C                      
  299. C     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX HERMITIAN             
  300. C     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING 
  301. C     REAL SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY  HTRIDI.                  
  302. C                      
  303. C     ON INPUT.        
  304. C                      
  305. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL                 
  306. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM                  
  307. C          DIMENSION STATEMENT.           
  308. C                      
  309. C        N IS THE ORDER OF THE MATRIX.    
  310. C                      
  311. C        AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY TRANS-                 
  312. C          FORMATIONS USED IN THE REDUCTION BY  HTRIDI  IN THEIR                
  313. C          FULL LOWER TRIANGLES EXCEPT FOR THE DIAGONAL OF AR.                  
  314. C                      
  315. C        TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS.            
  316. C                      
  317. C        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED.                
  318. C                      
  319. C        ZR CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED 
  320. C          IN ITS FIRST M COLUMNS.        
  321. C                      
  322. C     ON OUTPUT.       
  323. C                      
  324. C        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,     
  325. C          RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS     
  326. C          IN THEIR FIRST M COLUMNS.      
  327. C                      
  328. C     NOTE THAT THE LAST COMPONENT OF EACH RETURNED VECTOR   
  329. C     IS REAL AND THAT VECTOR EUCLIDEAN NORMS ARE PRESERVED. 
  330. C                      
  331. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,                
  332. C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY                 
  333. C                      
  334. C     ------------------------------------------------------------------        
  335. C                      
  336.       IF (M .EQ. 0) GO TO 200             
  337. C     .......... TRANSFORM THE EIGENVECTORS OF THE REAL SYMMETRIC               
  338. C                TRIDIAGONAL MATRIX TO THOSE OF THE HERMITIAN                   
  339. C                TRIDIAGONAL MATRIX. ..........              
  340.       DO 50 K = 1, N   
  341. C                      
  342.          DO 50 J = 1, M                   
  343.             ZI(K,J) = FLOP(-ZR(K,J)*TAU(2,K))                
  344.             ZR(K,J) = FLOP(ZR(K,J)*TAU(1,K))                 
  345.    50 CONTINUE         
  346. C                      
  347.       IF (N .EQ. 1) GO TO 200             
  348. C     .......... RECOVER AND APPLY THE HOUSEHOLDER MATRICES ..........          
  349.       DO 140 I = 2, N                     
  350.          L = I - 1     
  351.          H = AI(I,I)   
  352.          IF (H .EQ. 0.0D0) GO TO 140      
  353. C                      
  354.          DO 130 J = 1, M                  
  355.             S = 0.0D0                     
  356.             SI = 0.0D0                    
  357. C                      
  358.             DO 110 K = 1, L               
  359.                S = FLOP(S + AR(I,K)*ZR(K,J) - AI(I,K)*ZI(K,J))                  
  360.                SI = FLOP(SI + AR(I,K)*ZI(K,J) + AI(I,K)*ZR(K,J))                
  361.   110       CONTINUE   
  362. C     .......... DOUBLE DIVISIONS AVOID POSSIBLE UNDERFLOW ..........           
  363.             S = FLOP((S/H)/H)             
  364.             SI = FLOP((SI/H)/H)           
  365. C                      
  366.             DO 120 K = 1, L               
  367.                ZR(K,J) = FLOP(ZR(K,J) - S*AR(I,K) - SI*AI(I,K))                 
  368.                ZI(K,J) = FLOP(ZI(K,J) - SI*AR(I,K) + S*AI(I,K))                 
  369.   120       CONTINUE   
  370. C                      
  371.   130    CONTINUE      
  372. C                      
  373.   140 CONTINUE         
  374. C                      
  375.   200 RETURN           
  376.       END              
  377.       SUBROUTINE IMTQL2(NM,N,D,E,Z,IERR,JOB)                 
  378. C                      
  379.       INTEGER I,J,K,L,M,N,II,NM,MML,IERR  
  380.       DOUBLE PRECISION D(N),E(N),Z(NM,N)  
  381.       DOUBLE PRECISION B,C,F,G,P,R,S      
  382.       DOUBLE PRECISION FLOP               
  383. C                      
  384. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE IMTQL2,           
  385. C     NUM. MATH. 12, 377-383(1968) BY MARTIN AND WILKINSON,  
  386. C     AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE.   
  387. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971).           
  388. C                      
  389. C     THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS 
  390. C     OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE IMPLICIT QL METHOD.              
  391. C     THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO   
  392. C     BE FOUND IF  TRED2  HAS BEEN USED TO REDUCE THIS       
  393. C     FULL MATRIX TO TRIDIAGONAL FORM.    
  394. C                      
  395. C     ON INPUT.        
  396. C                      
  397. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL                 
  398. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM                  
  399. C          DIMENSION STATEMENT.           
  400. C                      
  401. C        N IS THE ORDER OF THE MATRIX.    
  402. C                      
  403. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.                  
  404. C                      
  405. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX                
  406. C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.    
  407. C                      
  408. C        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE                   
  409. C          REDUCTION BY  TRED2, IF PERFORMED.  IF THE EIGENVECTORS              
  410. C          OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN                
  411. C          THE IDENTITY MATRIX.           
  412. C                      
  413. C      ON OUTPUT.      
  414. C                      
  415. C        D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN                  
  416. C          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT                  
  417. C          UNORDERED FOR INDICES 1,2,...,IERR-1.             
  418. C                      
  419. C        E HAS BEEN DESTROYED.            
  420. C                      
  421. C        Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC                   
  422. C          TRIDIAGONAL (OR FULL) MATRIX.  IF AN ERROR EXIT IS MADE,             
  423. C          Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED               
  424. C          EIGENVALUES.                   
  425. C                      
  426. C        IERR IS SET TO                   
  427. C          ZERO       FOR NORMAL RETURN,  
  428. C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN    
  429. C  DETERMINED AFTER 30 ITERATIONS.        
  430. C                      
  431. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,                
  432. C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY                 
  433. C                      
  434. C     ------------------------------------------------------------------        
  435. C                      
  436. C                      
  437. C*****                 
  438. C     MODIFIED BY C. MOLER TO ELIMINATE MACHEP 11/22/78      
  439. C     MODIFIED TO ADD JOB PARAMETER 08/27/79                 
  440. C*****                 
  441.       IERR = 0         
  442.       IF (N .EQ. 1) GO TO 1001            
  443. C                      
  444.       DO 100 I = 2, N                     
  445.   100 E(I-1) = E(I)    
  446. C                      
  447.       E(N) = 0.0D0     
  448. C                      
  449.       DO 240 L = 1, N                     
  450.          J = 0         
  451. C     .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT ..........                 
  452.   105    DO 110 M = L, N                  
  453.             IF (M .EQ. N) GO TO 120       
  454. C*****                 
  455.             P = FLOP(DABS(D(M)) + DABS(D(M+1)))              
  456.             S = FLOP(P + DABS(E(M)))      
  457.             IF (P .EQ. S) GO TO 120       
  458. C*****                 
  459.   110    CONTINUE      
  460. C                      
  461.   120    P = D(L)      
  462.          IF (M .EQ. L) GO TO 240          
  463.          IF (J .EQ. 30) GO TO 1000        
  464.          J = J + 1     
  465. C     .......... FORM SHIFT ..........    
  466.          G = FLOP((D(L+1) - P)/(2.0D0*E(L)))                 
  467.          R = FLOP(DSQRT(G*G+1.0D0))       
  468.          G = FLOP(D(M) - P + E(L)/(G + DSIGN(R,G)))          
  469.          S = 1.0D0     
  470.          C = 1.0D0     
  471.          P = 0.0D0     
  472.          MML = M - L   
  473. C     .......... FOR I=M-1 STEP -1 UNTIL L DO -- ..........  
  474.          DO 200 II = 1, MML               
  475.             I = M - II                    
  476.             F = FLOP(S*E(I))              
  477.             B = FLOP(C*E(I))              
  478.             IF (DABS(F) .LT. DABS(G)) GO TO 150              
  479.             C = FLOP(G/F)                 
  480.             R = FLOP(DSQRT(C*C+1.0D0))    
  481.             E(I+1) = FLOP(F*R)            
  482.             S = FLOP(1.0D0/R)             
  483.             C = FLOP(C*S)                 
  484.             GO TO 160                     
  485.   150       S = FLOP(F/G)                 
  486.             R = FLOP(DSQRT(S*S+1.0D0))    
  487.             E(I+1) = FLOP(G*R)            
  488.             C = FLOP(1.0D0/R)             
  489.             S = FLOP(S*C)                 
  490.   160       G = FLOP(D(I+1) - P)          
  491.             R = FLOP((D(I) - G)*S + 2.0D0*C*B)               
  492.             P = FLOP(S*R)                 
  493.             D(I+1) = G + P                
  494.             G = FLOP(C*R - B)             
  495.             IF (JOB .EQ. 0) GO TO 185     
  496. C     .......... FORM VECTOR ..........   
  497.             DO 180 K = 1, N               
  498.                F = Z(K,I+1)               
  499.                Z(K,I+1) = FLOP(S*Z(K,I) + C*F)               
  500.                Z(K,I) = FLOP(C*Z(K,I) - S*F)                 
  501.   180       CONTINUE   
  502.   185       CONTINUE   
  503. C                      
  504.   200    CONTINUE      
  505. C                      
  506.          D(L) = FLOP(D(L) - P)            
  507.          E(L) = G      
  508.          E(M) = 0.0D0                     
  509.          GO TO 105     
  510.   240 CONTINUE         
  511. C     .......... ORDER EIGENVALUES AND EIGENVECTORS ..........                  
  512.       DO 300 II = 2, N                    
  513.          I = II - 1    
  514.          K = I         
  515.          P = D(I)      
  516. C                      
  517.          DO 260 J = II, N                 
  518.             IF (D(J) .GE. P) GO TO 260    
  519.             K = J      
  520.             P = D(J)   
  521.   260    CONTINUE      
  522. C                      
  523.          IF (K .EQ. I) GO TO 300          
  524.          D(K) = D(I)   
  525.          D(I) = P      
  526. C                      
  527.          IF (JOB .EQ. 0) GO TO 285        
  528.          DO 280 J = 1, N                  
  529.             P = Z(J,I)                    
  530.             Z(J,I) = Z(J,K)               
  531.             Z(J,K) = P                    
  532.   280    CONTINUE      
  533.   285    CONTINUE      
  534. C                      
  535.   300 CONTINUE         
  536. C                      
  537.       GO TO 1001       
  538. C     .......... SET ERROR -- NO CONVERGENCE TO AN           
  539. C                EIGENVALUE AFTER 30 ITERATIONS ..........   
  540.  1000 IERR = L         
  541.  1001 RETURN           
  542.       END              
  543.       SUBROUTINE CORTH(NM,N,LOW,IGH,AR,AI,ORTR,ORTI)         
  544. C                      
  545.       INTEGER I,J,M,N,II,JJ,LA,MP,NM,IGH,KP1,LOW             
  546.       DOUBLE PRECISION AR(NM,N),AI(NM,N),ORTR(IGH),ORTI(IGH) 
  547.       DOUBLE PRECISION F,G,H,FI,FR,SCALE  
  548.       DOUBLE PRECISION FLOP,PYTHAG        
  549. C                      
  550. C     THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF                 
  551. C     THE ALGOL PROCEDURE ORTHES, NUM. MATH. 12, 349-368(1968)                  
  552. C     BY MARTIN AND WILKINSON.            
  553. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).           
  554. C                      
  555. C     GIVEN A COMPLEX GENERAL MATRIX, THIS SUBROUTINE        
  556. C     REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS       
  557. C     LOW THROUGH IGH TO UPPER HESSENBERG FORM BY            
  558. C     UNITARY SIMILARITY TRANSFORMATIONS.                    
  559. C                      
  560. C     ON INPUT.        
  561. C                      
  562. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL                 
  563. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM                  
  564. C          DIMENSION STATEMENT.           
  565. C                      
  566. C        N IS THE ORDER OF THE MATRIX.    
  567. C                      
  568. C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING                   
  569. C          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED,   
  570. C          SET LOW=1, IGH=N.              
  571. C                      
  572. C        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,     
  573. C          RESPECTIVELY, OF THE COMPLEX INPUT MATRIX.        
  574. C                      
  575. C     ON OUTPUT.       
  576. C                      
  577. C        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,     
  578. C          RESPECTIVELY, OF THE HESSENBERG MATRIX.  INFORMATION                 
  579. C          ABOUT THE UNITARY TRANSFORMATIONS USED IN THE REDUCTION              
  580. C          IS STORED IN THE REMAINING TRIANGLES UNDER THE    
  581. C          HESSENBERG MATRIX.             
  582. C                      
  583. C        ORTR AND ORTI CONTAIN FURTHER INFORMATION ABOUT THE 
  584. C          TRANSFORMATIONS.  ONLY ELEMENTS LOW THROUGH IGH ARE USED.            
  585. C                      
  586. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,                
  587. C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY                 
  588. C                      
  589. C     ------------------------------------------------------------------        
  590. C                      
  591.       LA = IGH - 1     
  592.       KP1 = LOW + 1    
  593.       IF (LA .LT. KP1) GO TO 200          
  594. C                      
  595.       DO 180 M = KP1, LA                  
  596.          H = 0.0D0     
  597.          ORTR(M) = 0.0D0                  
  598.          ORTI(M) = 0.0D0                  
  599.          SCALE = 0.0D0                    
  600. C     .......... SCALE COLUMN (ALGOL TOL THEN NOT NEEDED) ..........            
  601.          DO 90 I = M, IGH                 
  602.    90    SCALE = FLOP(SCALE + DABS(AR(I,M-1)) + DABS(AI(I,M-1)))                
  603. C                      
  604.          IF (SCALE .EQ. 0.0D0) GO TO 180  
  605.          MP = M + IGH                     
  606. C     .......... FOR I=IGH STEP -1 UNTIL M DO -- ..........  
  607.          DO 100 II = M, IGH               
  608.             I = MP - II                   
  609.             ORTR(I) = FLOP(AR(I,M-1)/SCALE)                  
  610.             ORTI(I) = FLOP(AI(I,M-1)/SCALE)                  
  611.             H = FLOP(H + ORTR(I)*ORTR(I) + ORTI(I)*ORTI(I))  
  612.   100    CONTINUE      
  613. C                      
  614.          G = FLOP(DSQRT(H))               
  615.          F = PYTHAG(ORTR(M),ORTI(M))      
  616.          IF (F .EQ. 0.0D0) GO TO 103      
  617.          H = FLOP(H + F*G)                
  618.          G = FLOP(G/F)                    
  619.          ORTR(M) = FLOP((1.0D0 + G)*ORTR(M))                 
  620.          ORTI(M) = FLOP((1.0D0 + G)*ORTI(M))                 
  621.          GO TO 105     
  622. C                      
  623.   103    ORTR(M) = G   
  624.          AR(M,M-1) = SCALE                
  625. C     .......... FORM (I-(U*UT)/H)*A ..........              
  626.   105    DO 130 J = M, N                  
  627.             FR = 0.0D0                    
  628.             FI = 0.0D0                    
  629. C     .......... FOR I=IGH STEP -1 UNTIL M DO -- ..........  
  630.             DO 110 II = M, IGH            
  631.                I = MP - II                
  632.                FR = FLOP(FR + ORTR(I)*AR(I,J) + ORTI(I)*AI(I,J))                
  633.                FI = FLOP(FI + ORTR(I)*AI(I,J) - ORTI(I)*AR(I,J))                
  634.   110       CONTINUE   
  635. C                      
  636.             FR = FLOP(FR/H)               
  637.             FI = FLOP(FI/H)               
  638. C                      
  639.             DO 120 I = M, IGH             
  640.                AR(I,J) = FLOP(AR(I,J) - FR*ORTR(I) + FI*ORTI(I))                
  641.                AI(I,J) = FLOP(AI(I,J) - FR*ORTI(I) - FI*ORTR(I))                
  642.   120       CONTINUE   
  643. C                      
  644.   130    CONTINUE      
  645. C     .......... FORM (I-(U*UT)/H)*A*(I-(U*UT)/H) .......... 
  646.          DO 160 I = 1, IGH                
  647.             FR = 0.0D0                    
  648.             FI = 0.0D0                    
  649. C     .......... FOR J=IGH STEP -1 UNTIL M DO -- ..........  
  650.             DO 140 JJ = M, IGH            
  651.                J = MP - JJ                
  652.                FR = FLOP(FR + ORTR(J)*AR(I,J) - ORTI(J)*AI(I,J))                
  653.                FI = FLOP(FI + ORTR(J)*AI(I,J) + ORTI(J)*AR(I,J))                
  654.   140       CONTINUE   
  655. C                      
  656.             FR = FLOP(FR/H)               
  657.             FI = FLOP(FI/H)               
  658. C                      
  659.             DO 150 J = M, IGH             
  660.                AR(I,J) = FLOP(AR(I,J) - FR*ORTR(J) - FI*ORTI(J))                
  661.                AI(I,J) = FLOP(AI(I,J) + FR*ORTI(J) - FI*ORTR(J))                
  662.   150       CONTINUE   
  663. C                      
  664.   160    CONTINUE      
  665. C                      
  666.          ORTR(M) = FLOP(SCALE*ORTR(M))    
  667.          ORTI(M) = FLOP(SCALE*ORTI(M))    
  668.          AR(M,M-1) = FLOP(-G*AR(M,M-1))   
  669.          AI(M,M-1) = FLOP(-G*AI(M,M-1))   
  670.   180 CONTINUE         
  671. C                      
  672.   200 RETURN           
  673.       END              
  674.       SUBROUTINE COMQR3(NM,N,LOW,IGH,ORTR,ORTI,HR,HI,WR,WI,ZR,ZI,IERR           
  675.      *                 ,JOB)              
  676. C*****                 
  677. C     MODIFICATION OF EISPACK COMQR2 TO ADD JOB PARAMETER    
  678. C     JOB = 0  OUTPUT H = SCHUR TRIANGULAR FORM, Z NOT USED  
  679. C         = 1  OUTPUT H = SCHUR FORM, Z = UNITARY SIMILARITY 
  680. C         = 2  SAME AS COMQR2             
  681. C         = 3  OUTPUT H = HESSENBERG FORM, Z = UNITARY SIMILARITY               
  682. C     ALSO ELIMINATE MACHEP               
  683. C     C. MOLER, 11/22/78 AND 09/14/80     
  684. C     OVERFLOW CONTROL IN EIGENVECTOR BACKSUBSTITUTION, 3/16/82                 
  685. C*****                 
  686. C                      
  687.       INTEGER I,J,K,L,M,N,EN,II,JJ,LL,NM,NN,IGH,IP1,         
  688.      X        ITN,ITS,LOW,LP1,ENM1,IEND,IERR                 
  689.       DOUBLE PRECISION HR(NM,N),HI(NM,N),WR(N),WI(N),ZR(NM,N),ZI(NM,N),         
  690.      X       ORTR(IGH),ORTI(IGH)          
  691.       DOUBLE PRECISION SI,SR,TI,TR,XI,XR,YI,YR,ZZI,ZZR,NORM  
  692.       DOUBLE PRECISION FLOP,PYTHAG        
  693. C                      
  694. C     THIS SUBROUTINE IS A TRANSLATION OF A UNITARY ANALOGUE OF THE             
  695. C     ALGOL PROCEDURE  COMLR2, NUM. MATH. 16, 181-204(1970) BY PETERS           
  696. C     AND WILKINSON.   
  697. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).           
  698. C     THE UNITARY ANALOGUE SUBSTITUTES THE QR ALGORITHM OF FRANCIS              
  699. C     (COMP. JOUR. 4, 332-345(1962)) FOR THE LR ALGORITHM.   
  700. C                      
  701. C     THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS 
  702. C     OF A COMPLEX UPPER HESSENBERG MATRIX BY THE QR         
  703. C     METHOD.  THE EIGENVECTORS OF A COMPLEX GENERAL MATRIX  
  704. C     CAN ALSO BE FOUND IF  CORTH  HAS BEEN USED TO REDUCE   
  705. C     THIS GENERAL MATRIX TO HESSENBERG FORM.                
  706. C                      
  707. C     ON INPUT.        
  708. C                      
  709. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL                 
  710. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM                  
  711. C          DIMENSION STATEMENT.           
  712. C                      
  713. C        N IS THE ORDER OF THE MATRIX.    
  714. C                      
  715. C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING                   
  716. C          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED,   
  717. C          SET LOW=1, IGH=N.              
  718. C                      
  719. C        ORTR AND ORTI CONTAIN INFORMATION ABOUT THE UNITARY TRANS-             
  720. C          FORMATIONS USED IN THE REDUCTION BY  CORTH, IF PERFORMED.            
  721. C          ONLY ELEMENTS LOW THROUGH IGH ARE USED.  IF THE EIGENVECTORS         
  722. C          OF THE HESSENBERG MATRIX ARE DESIRED, SET ORTR(J) AND                
  723. C          ORTI(J) TO 0.0D0 FOR THESE ELEMENTS.              
  724. C                      
  725. C        HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS,     
  726. C          RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX.                
  727. C          THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN FURTHER          
  728. C          INFORMATION ABOUT THE TRANSFORMATIONS WHICH WERE USED IN THE         
  729. C          REDUCTION BY  CORTH, IF PERFORMED.  IF THE EIGENVECTORS OF           
  730. C          THE HESSENBERG MATRIX ARE DESIRED, THESE ELEMENTS MAY BE             
  731. C          ARBITRARY.                     
  732. C                      
  733. C     ON OUTPUT.       
  734. C                      
  735. C        ORTR, ORTI, AND THE UPPER HESSENBERG PORTIONS OF HR AND HI             
  736. C          HAVE BEEN DESTROYED.           
  737. C                      
  738. C        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,     
  739. C          RESPECTIVELY, OF THE EIGENVALUES.  IF AN ERROR    
  740. C          EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT   
  741. C          FOR INDICES IERR+1,...,N.      
  742. C                      
  743. C        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,     
  744. C          RESPECTIVELY, OF THE EIGENVECTORS.  THE EIGENVECTORS                 
  745. C          ARE UNNORMALIZED.  IF AN ERROR EXIT IS MADE, NONE OF                 
  746. C          THE EIGENVECTORS HAS BEEN FOUND.                  
  747. C                      
  748. C        IERR IS SET TO                   
  749. C          ZERO       FOR NORMAL RETURN,  
  750. C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN    
  751. C  DETERMINED AFTER A TOTAL OF 30*N ITERATIONS.              
  752. C                      
  753. C     MODIFIED TO GET RID OF ALL COMPLEX ARITHMETIC, C. MOLER, 6/27/79.         
  754. C                      
  755. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,                
  756. C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY                 
  757. C                      
  758. C     ------------------------------------------------------------------        
  759. C                      
  760.       IERR = 0         
  761. C*****                 
  762.       IF (JOB .EQ. 0) GO TO 150           
  763. C*****                 
  764. C     .......... INITIALIZE EIGENVECTOR MATRIX ..........    
  765.       DO 100 I = 1, N                     
  766. C                      
  767.          DO 100 J = 1, N                  
  768.             ZR(I,J) = 0.0D0               
  769.             ZI(I,J) = 0.0D0               
  770.             IF (I .EQ. J) ZR(I,J) = 1.0D0                    
  771.   100 CONTINUE         
  772. C     .......... FORM THE MATRIX OF ACCUMULATED TRANSFORMATIONS                 
  773. C                FROM THE INFORMATION LEFT BY CORTH ..........                  
  774.       IEND = IGH - LOW - 1                
  775.       IF (IEND) 180, 150, 105             
  776. C     .......... FOR I=IGH-1 STEP -1 UNTIL LOW+1 DO -- ..........               
  777.   105 DO 140 II = 1, IEND                 
  778.          I = IGH - II                     
  779.          IF (ORTR(I) .EQ. 0.0D0 .AND. ORTI(I) .EQ. 0.0D0) GO TO 140             
  780.          IF (HR(I,I-1) .EQ. 0.0D0 .AND. HI(I,I-1) .EQ. 0.0D0) GO TO 140         
  781. C     .......... NORM BELOW IS NEGATIVE OF H FORMED IN CORTH ..........         
  782.          NORM = FLOP(HR(I,I-1)*ORTR(I) + HI(I,I-1)*ORTI(I))  
  783.          IP1 = I + 1   
  784. C                      
  785.          DO 110 K = IP1, IGH              
  786.             ORTR(K) = HR(K,I-1)           
  787.             ORTI(K) = HI(K,I-1)           
  788.   110    CONTINUE      
  789. C                      
  790.          DO 130 J = I, IGH                
  791.             SR = 0.0D0                    
  792.             SI = 0.0D0                    
  793. C                      
  794.             DO 115 K = I, IGH             
  795.                SR = FLOP(SR + ORTR(K)*ZR(K,J) + ORTI(K)*ZI(K,J))                
  796.                SI = FLOP(SI + ORTR(K)*ZI(K,J) - ORTI(K)*ZR(K,J))                
  797.   115       CONTINUE   
  798. C                      
  799.             SR = FLOP(SR/NORM)            
  800.             SI = FLOP(SI/NORM)            
  801. C                      
  802.             DO 120 K = I, IGH             
  803.                ZR(K,J) = FLOP(ZR(K,J) + SR*ORTR(K) - SI*ORTI(K))                
  804.                ZI(K,J) = FLOP(ZI(K,J) + SR*ORTI(K) + SI*ORTR(K))                
  805.   120       CONTINUE   
  806. C                      
  807.   130    CONTINUE      
  808. C                      
  809.   140 CONTINUE         
  810. C*****                 
  811.       IF (JOB .EQ. 3) GO TO 1001          
  812. C*****                 
  813. C     .......... CREATE REAL SUBDIAGONAL ELEMENTS .......... 
  814.   150 L = LOW + 1      
  815. C                      
  816.       DO 170 I = L, IGH                   
  817.          LL = MIN0(I+1,IGH)               
  818.          IF (HI(I,I-1) .EQ. 0.0D0) GO TO 170                 
  819.          NORM = PYTHAG(HR(I,I-1),HI(I,I-1))                  
  820.          YR = FLOP(HR(I,I-1)/NORM)        
  821.          YI = FLOP(HI(I,I-1)/NORM)        
  822.          HR(I,I-1) = NORM                 
  823.          HI(I,I-1) = 0.0D0                
  824. C                      
  825.          DO 155 J = I, N                  
  826.             SI = FLOP(YR*HI(I,J) - YI*HR(I,J))               
  827.             HR(I,J) = FLOP(YR*HR(I,J) + YI*HI(I,J))          
  828.             HI(I,J) = SI                  
  829.   155    CONTINUE      
  830. C                      
  831.          DO 160 J = 1, LL                 
  832.             SI = FLOP(YR*HI(J,I) + YI*HR(J,I))               
  833.             HR(J,I) = FLOP(YR*HR(J,I) - YI*HI(J,I))          
  834.             HI(J,I) = SI                  
  835.   160    CONTINUE      
  836. C*****                 
  837.          IF (JOB .EQ. 0) GO TO 170        
  838. C*****                 
  839.          DO 165 J = LOW, IGH              
  840.             SI = FLOP(YR*ZI(J,I) + YI*ZR(J,I))               
  841.             ZR(J,I) = FLOP(YR*ZR(J,I) - YI*ZI(J,I))          
  842.             ZI(J,I) = SI                  
  843.   165    CONTINUE      
  844. C                      
  845.   170 CONTINUE         
  846. C     .......... STORE ROOTS ISOLATED BY CBAL ..........     
  847.   180 DO 200 I = 1, N                     
  848.          IF (I .GE. LOW .AND. I .LE. IGH) GO TO 200          
  849.          WR(I) = HR(I,I)                  
  850.          WI(I) = HI(I,I)                  
  851.   200 CONTINUE         
  852. C                      
  853.       EN = IGH         
  854.       TR = 0.0D0       
  855.       TI = 0.0D0       
  856.       ITN = 30*N       
  857. C     .......... SEARCH FOR NEXT EIGENVALUE ..........       
  858.   220 IF (EN .LT. LOW) GO TO 680          
  859.       ITS = 0          
  860.       ENM1 = EN - 1    
  861. C     .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT  
  862. C                FOR L=EN STEP -1 UNTIL LOW DO -- .......... 
  863.   240 DO 260 LL = LOW, EN                 
  864.          L = EN + LOW - LL                
  865.          IF (L .EQ. LOW) GO TO 300        
  866. C*****                 
  867.          XR = FLOP(DABS(HR(L-1,L-1)) + DABS(HI(L-1,L-1))     
  868.      X             + DABS(HR(L,L)) +DABS(HI(L,L)))           
  869.          YR = FLOP(XR + DABS(HR(L,L-1)))  
  870.          IF (XR .EQ. YR) GO TO 300        
  871. C*****                 
  872.   260 CONTINUE         
  873. C     .......... FORM SHIFT ..........    
  874.   300 IF (L .EQ. EN) GO TO 660            
  875.       IF (ITN .EQ. 0) GO TO 1000          
  876.       IF (ITS .EQ. 10 .OR. ITS .EQ. 20) GO TO 320            
  877.       SR = HR(EN,EN)   
  878.       SI = HI(EN,EN)   
  879.       XR = FLOP(HR(ENM1,EN)*HR(EN,ENM1))  
  880.       XI = FLOP(HI(ENM1,EN)*HR(EN,ENM1))  
  881.       IF (XR .EQ. 0.0D0 .AND. XI .EQ. 0.0D0) GO TO 340       
  882.       YR = FLOP((HR(ENM1,ENM1) - SR)/2.0D0)                  
  883.       YI = FLOP((HI(ENM1,ENM1) - SI)/2.0D0)                  
  884.       CALL WSQRT(YR**2-YI**2+XR,2.0D0*YR*YI+XI,ZZR,ZZI)      
  885.       IF (YR*ZZR + YI*ZZI .GE. 0.0D0) GO TO 310              
  886.       ZZR = -ZZR       
  887.       ZZI = -ZZI       
  888.   310 CALL WDIV(XR,XI,YR+ZZR,YI+ZZI,ZZR,ZZI)                 
  889.       SR = FLOP(SR - ZZR)                 
  890.       SI = FLOP(SI - ZZI)                 
  891.       GO TO 340        
  892. C     .......... FORM EXCEPTIONAL SHIFT ..........           
  893.   320 SR = FLOP(DABS(HR(EN,ENM1)) + DABS(HR(ENM1,EN-2)))     
  894.       SI = 0.0D0       
  895. C                      
  896.   340 DO 360 I = LOW, EN                  
  897.          HR(I,I) = FLOP(HR(I,I) - SR)     
  898.          HI(I,I) = FLOP(HI(I,I) - SI)     
  899.   360 CONTINUE         
  900. C                      
  901.       TR = FLOP(TR + SR)                  
  902.       TI = FLOP(TI + SI)                  
  903.       ITS = ITS + 1    
  904.       ITN = ITN - 1    
  905. C     .......... REDUCE TO TRIANGLE (ROWS) ..........        
  906.       LP1 = L + 1      
  907. C                      
  908.       DO 500 I = LP1, EN                  
  909.          SR = HR(I,I-1)                   
  910.          HR(I,I-1) = 0.0D0                
  911.          NORM = FLOP(DABS(HR(I-1,I-1)) + DABS(HI(I-1,I-1)) + DABS(SR))          
  912.          NORM = FLOP(NORM*DSQRT((HR(I-1,I-1)/NORM)**2 +      
  913.      X  (HI(I-1,I-1)/NORM)**2 + (SR/NORM)**2))               
  914.          XR = FLOP(HR(I-1,I-1)/NORM)      
  915.          WR(I-1) = XR                     
  916.          XI = FLOP(HI(I-1,I-1)/NORM)      
  917.          WI(I-1) = XI                     
  918.          HR(I-1,I-1) = NORM               
  919.          HI(I-1,I-1) = 0.0D0              
  920.          HI(I,I-1) = FLOP(SR/NORM)        
  921. C                      
  922.          DO 490 J = I, N                  
  923.             YR = HR(I-1,J)                
  924.             YI = HI(I-1,J)                
  925.             ZZR = HR(I,J)                 
  926.             ZZI = HI(I,J)                 
  927.             HR(I-1,J) = FLOP(XR*YR + XI*YI + HI(I,I-1)*ZZR)  
  928.             HI(I-1,J) = FLOP(XR*YI - XI*YR + HI(I,I-1)*ZZI)  
  929.             HR(I,J) = FLOP(XR*ZZR - XI*ZZI - HI(I,I-1)*YR)   
  930.             HI(I,J) = FLOP(XR*ZZI + XI*ZZR - HI(I,I-1)*YI)   
  931.   490    CONTINUE      
  932. C                      
  933.   500 CONTINUE         
  934. C                      
  935.       SI = HI(EN,EN)   
  936.       IF (SI .EQ. 0.0D0) GO TO 540        
  937.       NORM = PYTHAG(HR(EN,EN),SI)         
  938.       SR = FLOP(HR(EN,EN)/NORM)           
  939.       SI = FLOP(SI/NORM)                  
  940.       HR(EN,EN) = NORM                    
  941.       HI(EN,EN) = 0.0D0                   
  942.       IF (EN .EQ. N) GO TO 540            
  943.       IP1 = EN + 1     
  944. C                      
  945.       DO 520 J = IP1, N                   
  946.          YR = HR(EN,J)                    
  947.          YI = HI(EN,J)                    
  948.          HR(EN,J) = FLOP(SR*YR + SI*YI)   
  949.          HI(EN,J) = FLOP(SR*YI - SI*YR)   
  950.   520 CONTINUE         
  951. C     .......... INVERSE OPERATION (COLUMNS) ..........      
  952.   540 DO 600 J = LP1, EN                  
  953.          XR = WR(J-1)                     
  954.          XI = WI(J-1)                     
  955. C                      
  956.          DO 580 I = 1, J                  
  957.             YR = HR(I,J-1)                
  958.             YI = 0.0D0                    
  959.             ZZR = HR(I,J)                 
  960.             ZZI = HI(I,J)                 
  961.             IF (I .EQ. J) GO TO 560       
  962.             YI = HI(I,J-1)                
  963.             HI(I,J-1) = FLOP(XR*YI + XI*YR + HI(J,J-1)*ZZI)  
  964.   560       HR(I,J-1) = FLOP(XR*YR - XI*YI + HI(J,J-1)*ZZR)  
  965.             HR(I,J) = FLOP(XR*ZZR + XI*ZZI - HI(J,J-1)*YR)   
  966.             HI(I,J) = FLOP(XR*ZZI - XI*ZZR - HI(J,J-1)*YI)   
  967.   580    CONTINUE      
  968. C*****                 
  969.          IF (JOB .EQ. 0) GO TO 600        
  970. C*****                 
  971.          DO 590 I = LOW, IGH              
  972.             YR = ZR(I,J-1)                
  973.             YI = ZI(I,J-1)                
  974.             ZZR = ZR(I,J)                 
  975.             ZZI = ZI(I,J)                 
  976.             ZR(I,J-1) = FLOP(XR*YR - XI*YI + HI(J,J-1)*ZZR)  
  977.             ZI(I,J-1) = FLOP(XR*YI + XI*YR + HI(J,J-1)*ZZI)  
  978.             ZR(I,J) = FLOP(XR*ZZR + XI*ZZI - HI(J,J-1)*YR)   
  979.             ZI(I,J) = FLOP(XR*ZZI - XI*ZZR - HI(J,J-1)*YI)   
  980.   590    CONTINUE      
  981. C                      
  982.   600 CONTINUE         
  983. C                      
  984.       IF (SI .EQ. 0.0D0) GO TO 240        
  985. C                      
  986.       DO 630 I = 1, EN                    
  987.          YR = HR(I,EN)                    
  988.          YI = HI(I,EN)                    
  989.          HR(I,EN) = FLOP(SR*YR - SI*YI)   
  990.          HI(I,EN) = FLOP(SR*YI + SI*YR)   
  991.   630 CONTINUE         
  992. C*****                 
  993.       IF (JOB .EQ. 0) GO TO 240           
  994. C*****                 
  995.       DO 640 I = LOW, IGH                 
  996.          YR = ZR(I,EN)                    
  997.          YI = ZI(I,EN)                    
  998.          ZR(I,EN) = FLOP(SR*YR - SI*YI)   
  999.          ZI(I,EN) = FLOP(SR*YI + SI*YR)   
  1000.   640 CONTINUE         
  1001. C                      
  1002.       GO TO 240        
  1003. C     .......... A ROOT FOUND ..........  
  1004.   660 HR(EN,EN) = FLOP(HR(EN,EN) + TR)    
  1005.       WR(EN) = HR(EN,EN)                  
  1006.       HI(EN,EN) = FLOP(HI(EN,EN) + TI)    
  1007.       WI(EN) = HI(EN,EN)                  
  1008.       EN = ENM1        
  1009.       GO TO 220        
  1010. C     .......... ALL ROOTS FOUND.  BACKSUBSTITUTE TO FIND    
  1011. C                VECTORS OF UPPER TRIANGULAR FORM .......... 
  1012. C                      
  1013. C*****  THE FOLLOWING SECTION CHANGED FOR OVERFLOW CONTROL   
  1014. C       C. MOLER, 3/16/82                 
  1015. C                      
  1016.   680 IF (JOB .NE. 2) GO TO 1001          
  1017. C                      
  1018.       NORM = 0.0D0     
  1019.       DO 720 I = 1, N                     
  1020.          DO 720 J = I, N                  
  1021.             TR = FLOP(DABS(HR(I,J))) + FLOP(DABS(HI(I,J)))   
  1022.             IF (TR .GT. NORM) NORM = TR   
  1023.   720 CONTINUE         
  1024.       IF (N .EQ. 1 .OR. NORM .EQ. 0.0D0) GO TO 1001          
  1025. C     .......... FOR EN=N STEP -1 UNTIL 2 DO -- ..........   
  1026.       DO 800 NN = 2, N                    
  1027.          EN = N + 2 - NN                  
  1028.          XR = WR(EN)   
  1029.          XI = WI(EN)   
  1030.          HR(EN,EN) = 1.0D0                
  1031.          HI(EN,EN) = 0.0D0                
  1032.          ENM1 = EN - 1                    
  1033. C     .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- .......... 
  1034.          DO 780 II = 1, ENM1              
  1035.             I = EN - II                   
  1036.             ZZR = 0.0D0                   
  1037.             ZZI = 0.0D0                   
  1038.             IP1 = I + 1                   
  1039.             DO 740 J = IP1, EN            
  1040.                ZZR = FLOP(ZZR + HR(I,J)*HR(J,EN) - HI(I,J)*HI(J,EN))            
  1041.                ZZI = FLOP(ZZI + HR(I,J)*HI(J,EN) + HI(I,J)*HR(J,EN))            
  1042.   740       CONTINUE   
  1043.             YR = FLOP(XR - WR(I))         
  1044.             YI = FLOP(XI - WI(I))         
  1045.             IF (YR .NE. 0.0D0 .OR. YI .NE. 0.0D0) GO TO 765  
  1046.                YR = NORM                  
  1047.   760          YR = FLOP(YR/100.0D0)      
  1048.                YI = FLOP(NORM + YR)       
  1049.                IF (YI .NE. NORM) GO TO 760                   
  1050.                YI = 0.0D0                 
  1051.   765       CONTINUE   
  1052.             CALL WDIV(ZZR,ZZI,YR,YI,HR(I,EN),HI(I,EN))       
  1053.             TR = FLOP(DABS(HR(I,EN))) + FLOP(DABS(HI(I,EN))) 
  1054.             IF (TR .EQ. 0.0D0) GO TO 780  
  1055.             IF (TR + 1.0D0/TR .GT. TR) GO TO 780             
  1056.             DO 770 J = I, EN              
  1057.                HR(J,EN) = FLOP(HR(J,EN)/TR)                  
  1058.                HI(J,EN) = FLOP(HI(J,EN)/TR)                  
  1059.   770       CONTINUE   
  1060.   780    CONTINUE      
  1061. C                      
  1062.   800 CONTINUE         
  1063. C*****                 
  1064. C     .......... END BACKSUBSTITUTION ..........             
  1065.       ENM1 = N - 1     
  1066. C     .......... VECTORS OF ISOLATED ROOTS ..........        
  1067.       DO  840 I = 1, ENM1                 
  1068.          IF (I .GE. LOW .AND. I .LE. IGH) GO TO 840          
  1069.          IP1 = I + 1   
  1070. C                      
  1071.          DO 820 J = IP1, N                
  1072.             ZR(I,J) = HR(I,J)             
  1073.             ZI(I,J) = HI(I,J)             
  1074.   820    CONTINUE      
  1075. C                      
  1076.   840 CONTINUE         
  1077. C     .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE   
  1078. C                VECTORS OF ORIGINAL FULL MATRIX.            
  1079. C                FOR J=N STEP -1 UNTIL LOW+1 DO -- ..........                   
  1080.       DO 880 JJ = LOW, ENM1               
  1081.          J = N + LOW - JJ                 
  1082.          M = MIN0(J,IGH)                  
  1083. C                      
  1084.          DO 880 I = LOW, IGH              
  1085.             ZZR = 0.0D0                   
  1086.             ZZI = 0.0D0                   
  1087. C                      
  1088.             DO 860 K = LOW, M             
  1089.                ZZR = FLOP(ZZR + ZR(I,K)*HR(K,J) - ZI(I,K)*HI(K,J))              
  1090.                ZZI = FLOP(ZZI + ZR(I,K)*HI(K,J) + ZI(I,K)*HR(K,J))              
  1091.   860       CONTINUE   
  1092. C                      
  1093.             ZR(I,J) = ZZR                 
  1094.             ZI(I,J) = ZZI                 
  1095.   880 CONTINUE         
  1096. C                      
  1097.       GO TO 1001       
  1098. C     .......... SET ERROR -- NO CONVERGENCE TO AN           
  1099. C                EIGENVALUE AFTER 30 ITERATIONS ..........   
  1100.  1000 IERR = EN        
  1101.  1001 RETURN           
  1102.       END              
  1103.       SUBROUTINE WSVDC(XR,XI,LDX,N,P,SR,SI,ER,EI,UR,UI,LDU,VR,VI,LDV,           
  1104.      *                 WORKR,WORKI,JOB,INFO)                 
  1105.       INTEGER LDX,N,P,LDU,LDV,JOB,INFO    
  1106.       DOUBLE PRECISION XR(LDX,1),XI(LDX,1),SR(1),SI(1),ER(1),EI(1),             
  1107.      *                 UR(LDU,1),UI(LDU,1),VR(LDV,1),VI(LDV,1),                 
  1108.      *                 WORKR(1),WORKI(1)  
  1109. C                      
  1110. C                      
  1111. C     WSVDC IS A SUBROUTINE TO REDUCE A DOUBLE-COMPLEX NXP MATRIX X BY          
  1112. C     UNITARY TRANSFORMATIONS U AND V TO DIAGONAL FORM.  THE 
  1113. C     DIAGONAL ELEMENTS S(I) ARE THE SINGULAR VALUES OF X.  THE                 
  1114. C     COLUMNS OF U ARE THE CORRESPONDING LEFT SINGULAR VECTORS,                 
  1115. C     AND THE COLUMNS OF V THE RIGHT SINGULAR VECTORS.       
  1116. C                      
  1117. C     ON ENTRY         
  1118. C                      
  1119. C         X         DOUBLE-COMPLEX(LDX,P), WHERE LDX.GE.N.   
  1120. C                   X CONTAINS THE MATRIX WHOSE SINGULAR VALUE                  
  1121. C                   DECOMPOSITION IS TO BE COMPUTED.  X IS   
  1122. C                   DESTROYED BY WSVDC.   
  1123. C                      
  1124. C         LDX       INTEGER.              
  1125. C                   LDX IS THE LEADING DIMENSION OF THE ARRAY X.                
  1126. C                      
  1127. C         N         INTEGER.              
  1128. C                   N IS THE NUMBER OF COLUMNS OF THE MATRIX X.                 
  1129. C                      
  1130. C         P         INTEGER.              
  1131. C                   P IS THE NUMBER OF ROWS OF THE MATRIX X. 
  1132. C                      
  1133. C         LDU       INTEGER.              
  1134. C                   LDU IS THE LEADING DIMENSION OF THE ARRAY U                 
  1135. C                   (SEE BELOW).          
  1136. C                      
  1137. C         LDV       INTEGER.              
  1138. C                   LDV IS THE LEADING DIMENSION OF THE ARRAY V                 
  1139. C                   (SEE BELOW).          
  1140. C                      
  1141. C         WORK      DOUBLE-COMPLEX(N).    
  1142. C                   WORK IS A SCRATCH ARRAY.                 
  1143. C                      
  1144. C         JOB       INTEGER.              
  1145. C                   JOB CONTROLS THE COMPUTATION OF THE SINGULAR                
  1146. C                   VECTORS.  IT HAS THE DECIMAL EXPANSION AB                   
  1147. C                   WITH THE FOLLOWING MEANING               
  1148. C                      
  1149. C     A.EQ.0    DO NOT COMPUTE THE LEFT SINGULAR             
  1150. C               VECTORS.                  
  1151. C     A.EQ.1    RETURN THE N LEFT SINGULAR VECTORS           
  1152. C               IN U.  
  1153. C     A.GE.2    RETURNS THE FIRST MIN(N,P)                   
  1154. C               LEFT SINGULAR VECTORS IN U.                  
  1155. C     B.EQ.0    DO NOT COMPUTE THE RIGHT SINGULAR            
  1156. C               VECTORS.                  
  1157. C     B.EQ.1    RETURN THE RIGHT SINGULAR VECTORS            
  1158. C               IN V.  
  1159. C                      
  1160. C     ON RETURN        
  1161. C                      
  1162. C         S         DOUBLE-COMPLEX(MM), WHERE MM=MIN(N+1,P). 
  1163. C                   THE FIRST MIN(N,P) ENTRIES OF S CONTAIN THE                 
  1164. C                   SINGULAR VALUES OF X ARRANGED IN DESCENDING                 
  1165. C                   ORDER OF MAGNITUDE.   
  1166. C                      
  1167. C         E         DOUBLE-COMPLEX(P).    
  1168. C                   E ORDINARILY CONTAINS ZEROS.  HOWEVER SEE THE               
  1169. C                   DISCUSSION OF INFO FOR EXCEPTIONS.       
  1170. C                      
  1171. C         U         DOUBLE-COMPLEX(LDU,K), WHERE LDU.GE.N.   
  1172. C                   IF JOBA.EQ.1 THEN K.EQ.N,                
  1173. C                   IF JOBA.EQ.2 THEN K.EQ.MIN(N,P).         
  1174. C                   U CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS.            
  1175. C                   U IS NOT REFERENCED IF JOBA.EQ.0.  IF N.LE.P                
  1176. C                   OR IF JOBA.GT.2, THEN U MAY BE IDENTIFIED WITH X            
  1177. C                   IN THE SUBROUTINE CALL.                  
  1178. C                      
  1179. C         V         DOUBLE-COMPLEX(LDV,P), WHERE LDV.GE.P.   
  1180. C                   V CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS.            
  1181. C                   V IS NOT REFERENCED IF JOBB.EQ.0.  IF P.LE.N,               
  1182. C                   THEN V MAY BE IDENTIFIED WHTH X IN THE   
  1183. C                   SUBROUTINE CALL.      
  1184. C                      
  1185. C         INFO      INTEGER.              
  1186. C                   THE SINGULAR VALUES (AND THEIR CORRESPONDING                
  1187. C                   SINGULAR VECTORS) S(INFO+1),S(INFO+2),...,S(M)              
  1188. C                   ARE CORRECT (HERE M=MIN(N,P)).  THUS IF  
  1189. C                   INFO.EQ.0, ALL THE SINGULAR VALUES AND THEIR                
  1190. C                   VECTORS ARE CORRECT.  IN ANY EVENT, THE MATRIX              
  1191. C                   B = CTRANS(U)*X*V IS THE BIDIAGONAL MATRIX                  
  1192. C                   WITH THE ELEMENTS OF S ON ITS DIAGONAL AND THE              
  1193. C                   ELEMENTS OF E ON ITS SUPER-DIAGONAL (CTRANS(U)              
  1194. C                   IS THE CONJUGATE-TRANSPOSE OF U).  THUS THE                 
  1195. C                   SINGULAR VALUES OF X AND B ARE THE SAME. 
  1196. C                      
  1197. C     LINPACK. THIS VERSION DATED 07/03/79 .                 
  1198. C     G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.               
  1199. C                      
  1200. C     WSVDC USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS.    
  1201. C                      
  1202. C     BLAS WAXPY,PYTHAG,WDOTCR,WDOTCI,WSCAL,WSWAP,WNRM2,RROTG                   
  1203. C     FORTRAN DABS,DIMAG,DMAX1            
  1204. C     FORTRAN MAX0,MIN0,MOD,DSQRT         
  1205. C                      
  1206. C     INTERNAL VARIABLES                  
  1207. C                      
  1208.       INTEGER I,ITER,J,JOBU,K,KASE,KK,L,LL,LLS,LM1,LP1,LS,LU,M,MAXIT,           
  1209.      *        MM,MM1,MP1,NCT,NCTP1,NCU,NRT,NRTP1             
  1210.       DOUBLE PRECISION PYTHAG,WDOTCR,WDOTCI,TR,TI,RR,RI      
  1211.       DOUBLE PRECISION B,C,CS,EL,EMM1,F,G,WNRM2,SCALE,SHIFT,SL,SM,SN,           
  1212.      *                 SMM1,T1,TEST,ZTEST,SMALL,FLOP         
  1213.       LOGICAL WANTU,WANTV                 
  1214. C                      
  1215.       DOUBLE PRECISION ZDUMR,ZDUMI        
  1216.       DOUBLE PRECISION CABS1              
  1217.       CABS1(ZDUMR,ZDUMI) = DABS(ZDUMR) + DABS(ZDUMI)         
  1218. C                      
  1219. C     SET THE MAXIMUM NUMBER OF ITERATIONS.                  
  1220. C                      
  1221.       MAXIT = 75       
  1222. C                      
  1223. C     SMALL NUMBER, ROUGHLY MACHINE EPSILON, USED TO AVOID UNDERFLOW            
  1224. C                      
  1225. SHAR_EOF
  1226. #    End of shell archive
  1227. exit 0
  1228. -- 
  1229. Bob Page, U of Lowell CS Dept.  page@swan.ulowell.edu  ulowell!page
  1230. Have five nice days.
  1231.